Lymph node metastasis‐dependent molecular classification in papillary thyroid carcinoma defines aggressive metastatic outgrowth

Dear Editor, Lymph node metastasis (LNM) is the most important prognostic factor and a crucial indicator in the development of treatment strategies for patients with papillary thyroid carcinoma (PTC).1,2 A vast majority (25%–60%) of patients with PTC undergo thyroidectomy with neck node dissection.3 However, researchers continue to debate whether LNM is leading to over-treatment of patients with PTC.1,4 Interestingly, we have noted that upon stratifying LNM risk in PTC and clinical outcomes according to distinct molecular characteristics, LNM showing enrichment of genes related to inflammation and epithelial– mesenchymal transition (EMT)weremore aggressive than metabolically adapted LNM, and this was applicable in other tumours. A flowchart of this study is displayed in Figure S1. For discovery, we ran differential gene expression analysis based on the negative binomial distribution (DESEQ2) for 97 patients without LNM [LNM (−)] and 195 patients with LNM [LNM (+)] from our PTC patient cohort (Table S1). One hundred four differentially expressed genes (DEGs) were upregulated in LNM (+), whereas the other 140 DEGs were downregulated (Figure 1A). In LNM (+) specimens, DAVID software revealed significant clustering of up-regulated genes related with EMT pathways (Figure S2A–C). Subsequent analysis of DEGs using the Kmeans unsupervised clustering algorithm was performed to define marked molecular subtypes, wherein K = 2 showed the greatest classification outcome with a cophenetic coefficient value of.736 (Figure 1B, Figures S3 and S4). To develop a predictive scoring model, we employed LASSO regression and 18 gene signatures that effectively deciphered molecular subtypes with an area under curve value of .97 (Figure 1C–E, Figures S5A–C and S6A,B). We then characterised the molecular characteristics of differently defined clusters by evaluating hallmark


Lymph node metastasis-dependent molecular classification in papillary thyroid carcinoma defines aggressive metastatic outgrowth
Dear Editor, Lymph node metastasis (LNM) is the most important prognostic factor and a crucial indicator in the development of treatment strategies for patients with papillary thyroid carcinoma (PTC). 1,2 A vast majority (25%-60%) of patients with PTC undergo thyroidectomy with neck node dissection. 3 However, researchers continue to debate whether LNM is leading to over-treatment of patients with PTC. 1,4 Interestingly, we have noted that upon stratifying LNM risk in PTC and clinical outcomes according to distinct molecular characteristics, LNM showing enrichment of genes related to inflammation and epithelialmesenchymal transition (EMT) were more aggressive than metabolically adapted LNM, and this was applicable in other tumours.
A flowchart of this study is displayed in Figure S1. For discovery, we ran differential gene expression analysis based on the negative binomial distribution (DESEQ2) for 97 patients without LNM [LNM (−)] and 195 patients with LNM [LNM (+)] from our PTC patient cohort (Table S1). One hundred four differentially expressed genes (DEGs) were upregulated in LNM (+), whereas the other 140 DEGs were downregulated ( Figure 1A). In LNM (+) specimens, DAVID software revealed significant clustering of up-regulated genes related with EMT pathways ( Figure S2A-C). Subsequent analysis of DEGs using the Kmeans unsupervised clustering algorithm was performed to define marked molecular subtypes, wherein K = 2 showed the greatest classification outcome with a cophenetic coefficient value of.736 ( Figure 1B, Figures S3 and  S4). To develop a predictive scoring model, we employed LASSO regression and 18 gene signatures that effectively deciphered molecular subtypes with an area under curve value of .97 ( Figure 1C gene sets of a molecular signature database using single sample gene set enrichment analysis (ssGSEA). 5 PTC patients from Cluster 1 in both our cohort and TCGA shared up-regulated patterns of metabolism-related hallmark pathways. PTC patients from Cluster 2 in both our cohort and TCGA exhibited common enrichment of EMT and immune response ( Figure 2A, Figure S7). These differences in enrichment were not noted when patient samples were compared according to LNM status ( Figure S8A,B).
We noticed a dominant proportion of LNM (+) patients (133 [80.6%]) in Cluster 2 in our cohort. Remarkably, 80 N1b (cancer with lateral neck node metastasis) patient samples (48.5%) were allocated in Cluster 2, while 32 (25.2%) N1a (cancer with central neck node metastasis) and 30 (23.6%) N1b patient samples were in Cluster 1. The odds ratio to observe LNM in Cluster 2 was 4.098 ( Figure 2B). Validation with an odds ratio of 4.691 was observed for LNM in Cluster 2, which indicates a strong clinical association between Cluster 2 and LNM ( Figure S9). We then set out to determine the prognostic significance of our signature genes related with disease progression or recurrence by integrating clinical survival data from the TCGA database. In doing so, we noted that patients with LNM and high LASSO scores in Cluster 2 (upper two thirds of the designated cluster) displayed significant differences in their prognosis, compared to those with low LASSO scores (lower two-thirds of the designated cluster) among LNM patients in Cluster 1 ( Figure 2C). Interestingly, aggressive outcomes were also noticed in N1a patients from Cluster 2 ( Figure 2D). This is a crucial issue, as it is not clear whether pathological N1a metastasis can serve as a risk factor in patients with PTC. 6 Furthermore, we evaluated thyroid differentiation scores and other important PTC-related oncogenic pathways and found corroborating evidence of biological aggressiveness in Cluster 2 ( Figure Table S2). Meanwhile, progression-free survival according to LNM status or the presence of BRAF mutation, a unique genomic marker of PTC, did not show prognostic significance ( Figure S10A,B).
Additionally, we compared demographic and genomic information from PTC patients between different molecular clusters, as well as multivariate linear regression with LASSO score, to seek for any potential clinical relevance (Tables 1 and S3). Notably, BRAF V600E mutation, T stage, degree of LNM, and extrathyroidal extension (ETE) differed significantly between two clusters, and depth of ETE exhibited a significant proportionate relationship with LASSO score, suggesting a greater chance of severe metastasis in Cluster 2 PTC patients (Tables S4-S6).
Our gene expression program was extended to actual LNM samples of PTC by applying it to public databases, GSE60542 and GSE151179 ( Figure 3A, S11A-F). 7,8 Also, we challenged the applicability of our gene expression program in LNM of other tumour types and incorporated datasets of breast cancer and melanoma from TCGA, GSE56493, and GSE65904. 9,10 Likewise, molecular subtypes were deciphered using our proposed LASSO gene expression program, and identified molecular characteristics were validated in both primary tumour and LNM ( Figure 3A). Moreover, LASSO scores of paired primary tumour and LNM samples showed significant relevance, suggesting preservation of molecular profiles after disease progression ( Figure S12A-E). Analogous tendencies were observed in poorly differentiated thyroid carcinoma and anaplastic thyroid cancer patients in GSE76039 ( Figure  S13A,B).
To expand the program's clinical significance, we analysed treatment responses to radioiodine therapy among PTC patients in GSE151179. Notably, significant radioactive   iodine (RAI) uptake was found in Cluster 1 with markedly lower LASSO scores than patients without radioiodine uptake ( Figure 3B,C). In addition, LASSO scores were higher in patients who underwent RAI treatment (RAIT) ( Figure 3D) and in an intermediate-or high-dose group, compared to no treatment or low-dose group ( Figure S14A, S14B). LASSO scores also showed a clear tendency to increase as suppressed Tg increased among all patients who underwent RAIT ( Figure 3E) and patients who underwent high-dose RAIT ( Figure S14C). Finally, we observed concordance between the aggressive outcomes noticed in Cluster 2 primary tumours with LNM and LNM samples from other tumour types ( Figure 3F-K).
Herein, we report two potential distinct LNM molecular subtypes of PTC with different clinical outcomes in terms of LNM and biological features. Proposed signatures were conserved throughout disease progression and validated in other thyroid cancer data and different tumour types. An EMT/inflammation axis, which we suggest as an aggressive LNM mechanism, requires support from additional studies for potential application in clinical settings.

A C K N O W L E D G M E N T S
We thank Ji Young Kim, Hee Chang Yu, and Hoyoung Kim for their excellent technical support.

C O N F L I C T S O F I N T E R E S T S S TAT E M E N T
The authors declare that they have no competing interests.

F I G U R E 3 Validation of molecular subtypes in LNM samples and other cancer cohorts. (A)
Results of ssGSEA analysis for GSE60542 paired primary tumour and LNM samples and GSE151179 primary tumour and LNM samples were figured into a heatmap. Results of enrichment analysis for TCGA-BRCA, GSE56493, TCGA-SKCM, and GSE65904 were enrolled for comparison. Genes are listed in the order of significance evaluated in our discovery cohort. Fold change was calculated as a comparison of average enrichment scores in each molecular cluster. (B) Panels representing RAI uptake and mutational status of PTC samples in GSE151179 are depicted as molecular clusters. Fisher's exact test was performed to evaluate differences in the frequency of RAI uptake in both clusters. (C) LASSO scores were compared between patients in GSE151179 with different responses after RAIT. (D, E) LASSO scores were compared in our in-house database of patients according to whether or not RAIT was performed and serum thyroglobulin levels 1 year after RAIT. Mann-Whitney U-test was performed for statistical analysis. (F, H, I, J) Disease-specific survival was compared between molecular subtypes in either PTC samples with LNM or LNM. Log rank (Mantel-Cox) test was performed to prove significance of Kaplan-Meier survival curves. (G) Overall survival rate was assessed between molecular subtypes of GSE56493 LNM samples. (K) Metastatic free survival was compared according to molecular subtypes from GSE65904 LNM samples. *P < .05, **P < .01, ***P < .001. BRCA, breast cancer; LNM, lymph node metastasis; PT, primary tumour; RAI, radioactive iodine treatment; RAI-R, radioactive iodine treatment refractory; SKCM, skin cutaneous melanoma; Tg, thyroglobulin.